Discovery of Localized Regions of Excess 10-TeV Cosmic Rays 



A. A. Abdo/ B. Allcn,^ T. Aunc,^ 
R. W. Ellsworth,^ L. Flcysher,^ R. 
P. H. HiintcmcycrJ B. E. Koltcrman,^ 



D. Berlcy,4 E. Blaufuss,-^ S. Casanova,^ C. Chcii,^ B. L. Dingus/ 
Fleysher,^ M. M. Gonzalez, 1° J. A. Goodman ^ C. M. Hoffman/ 
C. P. LansdcU," J. T. Linncmann/^ J. E. McEnery/^ A. I. Mincer/ 



OO 
O 

o 

o 
O 



6 



> 

00 

rn 

O 

OO 

O 



X 



Ncmcthy/ D 
J. Smith/ G 



Noycs/ J. Prctz/ J. M. Ryan/"* P. M. Saz Parkinson/ A. Shoup/^ G. Sinnis/ 

W. Sullivan/ V. Vasilciou/ G. P. Walker/^ D. A. Williams/ and G. B. Yodh^ 

'Naval Research Laboratory, Washington, DC 
''Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 
'^University of California, Santa Cruz, CA 
^ University of Maryland, College Park, MD 
^Max Planck Institut fur Kernphysik, Heidelberg, Germany 
University of California, Irvine, CA 
'Los Alamos National Laboratory, Los Alamos, NM 
^George Mason University, Fairfax, VA 
^New York University, New York, NY 
'" Instituto de Astronomia, Universidad Nacional Autonoma de Mexico, D.F., MEXICO 
'^'^ Institute for Defense Analyses, Alexandria, VA 
'^'^ Michigan State University, East Lansing, MI 
'^NASA Goddard Space Flight Center, Greenbelt, MD 
University of New Hampshire, Durham, NH 
'^Ohio State University, Lima, OH 
National Security Technologies, Las Vegas, NV 

The 7 year data set of the Milagro TeV observatory contains 2.2 x 10^^ events of which most 
are due to hadronic cosmic rays. This data is searched for evidence of intermediate scale structure. 
Excess emission on angular scales of ~ 10° has been found in two localized regions of unknown origin 
with greater than 12a significance. Both regions are inconsistent with pure gamma-ray emission 
with high confidence. One of the regions has a difi'erent energy spectrum than the isotropic cosmic- 
ray flux at a level of 4.6(j, and it is consistent with hard spectrum protons with an exponential 
cutoff, with the most significant excess at ~ 10 TeV. Potential causes of these excesses are explored, 
but no compelling explanations are found. 

PACS numbers: 95.55. Vj, 95.85.Ry, 96.50.Xy, 98.35. Eg, 98.70.Sa 
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The flux of charged cosmic rays at TeV energies is 
known to be nearly isotropic. This is due to Galac- 
tic magnetic fields, which randomize the directions 
of charged particles. However, numerous experiments 
across a wide range of energies have found anisotropy on 
large angular scales, typically with a fractional amplitude 
of ~ 10""^ (see [H, 0, H, 0, m , for example). Large-scale 
anisotropy is also seen in data from the Milagro detector 
, here we present the results of an analysis sensitive to 
intermediate angular scales (~ 10°). 

Milagro Q is a water Cherenkov air shower detector 
located in New Mexico, USA at an altitude of 2630m 
and at 36° N latitude. It is composed of a central 60m x 
80m pond surrounded by a sparse 200m x 200m array of 
175 "outrigger" water tanks. The pond is instrumented 
with 723 photomultiplier tubes (PMTs) in two layers. 
The top layer and outrigger tanks are used to determine 
the direction and energy, while the bottom layer is used 
to distinguish between gamma-ray induced and hadron 
induced air showers. The outriggers, with each tank con- 
taining a single PMT, improve the angular and energy 
resolution of the detector for events collected after May, 



2003. Milagro has a ^2 sr fleld of view, operates with 
a >90% duty cycle, and has a trigger rate from cosmic 
rays of ~1700 Hz, making it well-suited to searching for 
anisotropy in the arrival directions of TeV cosmic rays. 

For studies on small to intermediate scales (< 10°), 
an adaptation of the gamma-ray point source analysis, 
which has been published previously 0, is used. The 
primary difference between the previous analysis and the 
current analysis is that no cosmic-ray background rejec- 
tion cuts are made. These cuts removed over 90% of the 
events, so the analysis reported here uses nearly 10 times 
the number of events of the previous analysis. Like the 
previous analysis, a signal map is made based on the ar- 
rival direction of each event. A background map is also 
created using the "direct integration" technique 0, in 
which two-hour intervals arc used to calculate the back- 
ground. Because of this two-hour interval, the analysis 
is relatively insensitive to features larger than ^ 30° in 
right ascension (RA); a different analysis of the Milagro 
data sensitive to larger features has been performed and 
is presented elsewhere 0. 

In the gamma-ray point source analysis, the signal and 
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FIG. 1: Map of significances for tlie Milagro data set without any cuts to remove tlie hadronic cosmic-ray background. A 
10° bin was used to smooth the data, and the color scale gives the significance. The solid line marks the Galactic plane, and 
every 10° in Galactic latitude are shown by the dashed lines. The black dot marks the direction of the heliotail, which is the 
direction opposite the motion of the solar system with respect to the local interstellar matter. The fractional excess of Region 
A is ~ 6 X 10~*, while for Region B it is ~ 4 x 10~*. The deep deficits bordering the regions of excess appear because the 
background calculation has been raised by the excess. 



background maps arc smoothed with a square bin of size 
2.1°/ cos((5) in RA by 2.1° in Declination (5), which is op- 
timal for Milagro's angular resolution. However, the bin 
size may be increased to improve the sensitivity to larger 
features, with a maximum size of about 10° for 6 < 60° 
(for 6 > 60°, the RA bin width 10°/cos(5) becomes too 
large for the 30° background interval). The significance 
is calculated using the method of Li and Ma [8| . 

The analysis has been applied to data collected be- 
tween July 2000 and August 2007. Events were required 
to have a zenith angle < 45° and nFit > 20, where 
nFit is the number of PMTs used in the angle fit. With 
these cuts, the dataset consists of 2.2 x 10^"'^ events with 
a median energy of 1 TeV and an average angular 
resolution of < 1°. Figure [1] shows the map of signif- 
icances made with 10° smoothing and no cuts to dis- 
criminate gamma rays from charged cosmic rays. The 
Cygnus Region, which has previously been shown to emit 
TeV gamma rays is clearly visible. The excesses la- 
beled "Region A" and "Region B" are seen with peak 
significances of 15. Oct and 12.7ct, respectively. These are 
pre-trial significances because the location and extent of 
the excesses were determined by examining the data. A 
map such as shown in Figured] has a few 100,000 inde- 
pendent bins, but given the high statistical significance 
many maps could be examined and the post trials signif- 
icance would be reduced by < Ict. The fractional excess 
relative to the cosmic-ray background is 6 x 10^"* for 
Region A and ~ 4 x 10~^ for Region B. Note that both 
excesses are paralleled by regions of deep deficit; this is a 
known effect of the analysis due to the fact that Regions 
A and B arc included in the background calculation of 
neighboring areas in RA. Therefore, the excess raises the 



background calculation above its actual value, resulting 
in an apparent deficit. 

Similarity is seen between the map in Figure [T] and 
results from the Tibet AS7 collaboration [1, H^Ij but a 
direct comparison cannot be made because the analysis 
methods differ. For each band in 5, the Tibet analysis 
measured the excess (or deficit) relative to the average 
for that 5 band, making it sensitive to the large-scale 
anisotropy discussed in Q. Smaller features, such as Re- 
gions A and B, were superimposed on the large-scale vari- 
ation, which is several times greater in amplitude. Con- 
versely, in the analysis presented here, the excess/deficit 
was measured with respect to the local background cal- 
culation, which is determined from the data ±30° in RA. 
This is illustrated in Figure [21 which shows the data and 
background calculation versus RA for a 10° band in dec- 
lination without any smoothing applied to the data. The 
large-scale variation dominates the figure, but the back- 
ground calculation makes the analysis sensitive only to 
features with an extent smaller than ~ 30° in RA. It is 
noteworthy that the Tibet AS 7 collaboration has devel- 
oped a model for the large-scale structure, and the resid- 
ual map after subtracting that model from their data 
shows excesses similar to Regions A and B [lo| . 

To estimate the extent of Region A, an elliptical Gaus- 
sian was fit to the excess map of the data in 0.1° bins 
prior to smoothing. The fit, which accounted for the 
change in sensitivity with declination, returned a cen- 
troid of RA = 69.4° ± 0.7°, 6 = 13.8° ± 0.7°, a half width 
of 2.6° ±0.3°, a half length of 7.6° ±1.1°, and an angle of 
46° ± 4° with respect to the RA axis. It is important to 
note that this fit focused on a "hot spot" in the general 
excess of Region A, but there is still excess extending to 
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FIG. 2: Signal and background events vs RA for 10° < 5 < 
20°. The plot was made using independent 10° 6 by 1° RA 
bins (i.e. no smoothing). A subset of the data was used 
in which there are only full days of data in order to give an 
approximately uniform exposure in RA. Region A corresponds 
to the excess at RA« 70°. This plot shows that the Region 
A excess is inherent in the raw signal data and is not due to 
an underestimation of the background. 

lower declinations. A fit was not performed to the excess 
in Region B due to its large, irregular shape. 

While the excesses in Regions A and B are statistically 
significant, systematic causes must be ruled out. Poten- 
tial weather-related effects were explored by dividing the 
data into the four seasons, and both excesses were seen 
in each season. The data were also divided into yearly 
datasets to investigate whether changes to the detector 
could play a role, and again the excesses were found in 
each datasct. The analysis was also run using universal 
time instead of sidereal time to check for day/night effects 
which could masquerade as a signal. In addition, the data 
were analyzed using anti-sidereal time, which provides a 
sanity check on the analysis since it will scramble real 
celestial features. No excess appears in either analysis. 

Potential errors in the background calculation were 
also investigated. Figure [2] shows the number of events 
versus RA for the signal and background for 10° < (5 < 
20°, using independent 10° 5 by 1° RA bins (i.e. no 
smoothing). The data for this figure were chosen to in- 
clude only full days in order to achieve an approximately 
uniform exposure as a function of RA, and the broad 
deficit seen by the Tibet Air Shower Array is evident 
(centered around RA = 180°). As can be seen, the back- 
ground estimate as calculated via the direct integration 
technique 0] agrees well with the data. The excess cor- 
responding to Region A is clearly inherent in the raw 
signal data and is not an artifact created by the back- 
ground subtraction. A similar result is found for Region 
B. 

Diagnostic tests have been performed to gain insight 
into the nature of Regions A and B. For the purposes 



of these tests, Region A is defined as the box bounded 
by 66° < RA < 76° and 10° < (5 < 20°. Region B is 
defined as the union of two boxes: 117° < RA < 131° and 
15° < S < 40°, and 131° < RA < 141° and 40° < 5 < 
50°. These definitions were chosen by visual inspection 
of Figure [TJ 

To check for flux variation, the analysis was applied 
to yearly and seasonal datasets. For each region, the 
yearly excess was consistent with a constant flux. Both 
regions also had a significant excess during each of the 
four seasons, with the respective fractional excess in parts 
per 10000 in spring, summer, fall, and winter of 6.5 ± 
0.9, 4.0 ± 0.9, 6.4 ± 0.9, 7.2 ± 0.9 for region A and 3.5 ± 
0.4, 3.3 ± 0.4, 4.0 ± 0.4, , 4.7 ± 0.4 for region B. In both 
cases the fractional excess was lowest in the summer and 
highest in the winter, and the probability relative to 
a constant fractional excess is only about 5% for each 
region. While this may provide insight into the cause of 
these excesses, only statistical errors are given. There 
could be systematic effects such as the slightly higher 
energy threshold of Milagro in winter when there is snow 
on top of the pond. 

The excesses in Regions A and B are inconsistent with 
pure gamma-ray emission. We can statistically separate 
gamma-ray events from cosmic-ray events utilizing two 
parameters. The compactness parameter^ uses PMT 
information in the bottom layer of Milagro to identify 
the penetrating particles characteristic of a hadronic air 
shower. The distribution of compactness depends on the 
energy spectrum of the source with higher energy gamma 
rays producing showers of greater compactness. In order 
to exclude a gamma-ray hypothesis of any spectrum, we 
also fit an energy parameter iout-, the fraction of outrigger 
PMTs that detect light. Figure [3^ shows the fractional 
excess of loge(foMt) for regions A and B. We hypothesize 
a spectrum for the excess of the form: 

dN/dE (xE^e-^- (1) 

where 7 is the spectral index and Ec is the characteris- 
tic energy at which the spectrum cuts off. We attempt 
this fit for regions A and B assuming separately that the 
primary particles are purely gamma rays and purely pro- 
tons. The gamma-ray hypothesis in region A has a 
of 124.0 with 16 degrees of freedom. The cosmic-ray hy- 
pothesis produces a reasonable fit, with a minimum of 
10.3. For region B, the best gamma-ray hypothesis has 
a of 84.8 compared to a best cosmic-ray hypothesis of 
~ 19.0, again with 16 degrees of freedom. Thus the 
proton hypothesis is a reasonable fit for both regions and 
the gamma-ray hypothesis is inconsistent with probabil- 
ities of 9 X 10^^^ for region A and 2 x 10^^^ for region B. 
The possibility that the regions contain some admixture 
of protons and gamma rays has not been considered. Fig- 
ure |4] shows the la, 2a, and ia regions around the best 
fit for region A and region B for the pure proton hypoth- 
esis. Some care must be taken in interpreting Figure 21 
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FIG. 3: (a): Differential plot of the fractional excess versus 
loge(fotit) for Regions A and B, where font is the fraction of the 
outriggers hit. The spectrum of Region A is significantly dif- 
ferent than the background (2 x 10^") , which is represented by 
the horizontal line, (b): Profile plot of the simulated energy 
of protons for the as a function of loge(foiit)- The ranges are 
asymmetric and contain the inner 68% of simulated events. 

It does not account for our systematic errors. There is 
a estimated systematic uncertainty in the spectral index 
of ±0.2 due to variation in the trigger threshold (caused 
by such things as changes in atmospheric pressure or ice 
on the pond). There is also a ~ 30% systematic un- 
certainty in the energy scale due to the threshold varia- 
tion, as well as discrepancy between the simulated and 
measured trigger rates. The fit does not constrain the 
spectrum well except to suggest that a hard spectrum 
is favored, particularly for region A. The cutoff energy 
is constrained, with logio{Ec/GeV) ~ 4-0-o5(stat) f'^'^ 
region A and logio(£'c/G'ey) = 4.0+°;5(^(^j) for region 
B. Most importantly, the pure gamma-ray hypothesis is 
strongly disfavored. 

We can see that the excesses in regions A and B are 
harder than the spectrum of the isotropic part of cosmic 
rays with minimal systematic effects by looking at the 
data alone. Figure [3^ shows the fractional excess in rc- 



FIG. 4: Results of a fit to the excesses in region A and 
B, assuming a pure-proton spectrum of the form in Equation 
[T] The top panel shows the results for region A and the bot- 
tom panel shows the results for region B. The la, 2a and So- 
allowed regions of the spectral index 7 and the cutoff energy 
Ec are indicated by the shaded regions. 

gions A and B as a function of font- Assuming that the 
excess is due to cosmic rays of the same spectrum, we 
would expect the fractional excess to be completely flat. 
The offset from zero tells us that this region is in fact 
an excess. A test of whether the curves in figure [3] 
are flat for region A(B) returns a chance probability of 
2 X 10~^ (6 X 10""^), independent of systematic errors. 
The excess of Region A is most significantly detected for 
loge(foMt) ^ —1.5, corresponding to an energy of about 
10 TeV for protons, as shown in Figure [3)d. At around 
10 TeV, the spectrum cuts off consistent with the results 
of the spectral fit. 

There is currently no compelling explanation for the 
excesses in Regions A and B. One possibility is that they 
could be due to neutrons, but this is unlikely because 
the decay length of 10 TeV neutrons is only about 0.1 
parsecs, which is much closer than the nearest star. An- 
other possibility is that these excesses could be caused 
by a Galactic cosmic-ray accelerator, but this is difficult 
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because the gyroradius of a 10 TeV proton in a 2/iG mag- 
netic field, which is the estimated strength of the local 



Galactic field [ll|, is only ^ 0.005 parsecs. In order for 



protons from a cosmic-ray accelerator to reach us, the 
intervening magnetic field must connect us to the source 
and be coherent out to ~ 100 parsecs since there are likely 
no sources within this distance. However, the direction 
of both regions is nearly perpendicular to the expected 
Galactic magnetic field direction [ll| With non-standard 
cosmic-ray diffusion, it is conceivable to account for these 
regions w^th a nearby cosmic-ray accelerators. 

Another possibility is that one or both of the excesses 
could be caused by the hcliospherc. This explanation is 
supported by the coincidence of Region A with the direc- 
tion of the hehotail (RA« 74°, 5 w 17° [i3|), which is the 
direction opposite the motion of the solar system with re- 
spect to the local interstellar matter. The possibility that 
we are seeing neutron production in the gravitationally- 
focused tail of inter-stellar medium material has been 
considered and discarded in [l2| because of insufficient 
target material. 

In summary, Milagro has observed two unexplained re- 
gions of excess with high significance. Potential system- 
atic causes have been examined and excluded. Both ex- 
cesses are inconsistent with pure gamma rays with high 
confidence, and their energy spectra arc moderately to 
strongly inconsistent with the spectrum of the isotropic 
cosmic-ray flux. In particular, the excess in Region A 
can be modeled as hard spectrum protons with a cutoff. 
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